# Country outliers
outliers = c('EGY', 'IRQ', 'AFG', 'ISR', 'COL', 'PAK')
f = oda ~ oda_lag + p_polity2 + gdppc + pop + factor(year)*factor(agency)
f1 = update(f, '.~.+alliance*independence')
f2 = update(f, '.~.+trade*independence')
f3 = update(f, '.~.+infant_mortality*independence')

mod1 = lm(f1, dat[!dat$iso3c %in% outliers,])
mod2 = lm(f2, dat[!dat$iso3c %in% outliers,])
mod3 = lm(f3, dat[!dat$iso3c %in% outliers,])

country_outliers = list(mod1, mod2, mod3)
se = lapply(country_outliers, function(i) sqrt(diag(sandwich(i))))

print(screenreg(country_outliers, omit.coef='factor', digits=3, override.se=se))

# Agency outliers
mod1 = lm(f1, dat[dat$agency!='AID',])
mod2 = lm(f2, dat[dat$agency!='AID',])
mod3 = lm(f3, dat[dat$agency!='AID',])
mod4 = lm(f1, dat[dat$agency!='ADF',])
mod5 = lm(f2, dat[dat$agency!='ADF',])
mod6 = lm(f3, dat[dat$agency!='ADF',])
mod7 = lm(f1, dat[dat$agency!='DOD',])
mod8 = lm(f2, dat[dat$agency!='DOD',])
mod9 = lm(f3, dat[dat$agency!='DOD',])

agency_outliers = list(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8, mod9)
se = lapply(agency_outliers, function(i) sqrt(diag(sandwich(i))))

print(screenreg(agency_outliers, omit.coef='factor', digits=3, override.se=se))
